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We investigate a simple variation of the Generalized Harmonic method for evolving the Einstein 
equations. A flat space wave equation for metric perturbations is separated from the Ricci tensor, 
with the rest of the Ricci tensor becoming a source for these wave equations. We demonstrate 
that this splitting method allows for the accurate simulation of compact objects, with gravitational 
field strengths less than or equal to those of neutron stars. This method could thus provide a 
straightforward path for general relativistic effects to be added to astrophysics simulations, such as 
in core collapse, accretion disks, and extreme mass ratio systems. 
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I. INTRODUCTION 

The wave equation is one of the best known and most 
studied partial differential equations, and it is solved by a 
variety of numerical methods in computer simulations in- 
volving radiative problems. Wave equations also appear 
as the leading order terms within the Einstein equations 
of General Relativity (GR), when selecting coordinates 
that follow the harmonic coordinate condition Dx"" — 0. 
The Generalized Harmonic (GH) formulation of GR ex- 
pands on the harmonic coordinate condition by setting 
the d'Alembertian of the coordinates to be a general func- 
tion: Dx"" = H°^. This preserves the Einstein equations 
in a wave equation form, but now allows for any coordi- 
nate system to be used. The GH formulation of GR was 
the first to enable the simulation of the inspiral, merger 
and ringdown of two black holes [l| (and wasquickly fol- 
lowed by groups using the BSSN formalism - for a 
recent comparison of the two methods see [4|). 

We report here on a simple rewriting of the GH for- 
mulation that further exploits the wave equation nature 
of the Einstein equations. Flat space wave equations for 
metric perturbations are split off from the Ricci tensor, 
and the rest of the Ricci tensor is converted into sources 
for these wave equations. The complexities of the Ein- 
stein field equations are thus packaged within these new 
source terms, which we find to be well behaved when the 
gravitational field strengths are on the order of those for 
a neutron star or weaker. Note that similar methods for 
splitting off flat space wave equations have existed for a 
long time, dating back to de Bonder: [ii (and see also 
0,1). 

To determine whether this splitting is useful numer- 
ically we perform a simple simulation of a binary neu- 
tron star inspiral and merger. We find that the gravi- 
tational fields evolve stably, that we can extract the ex- 
pected gravitational waves, and we verify that the Hamil- 
tonian constraint violations converge to zero. This test 
shows that this splitting method provides a fairly simple 
way for general relativistic effects and gravitational wave 



production to be added to astrophysics simulations that 
have traditionally only used Newtonian gravity or sim- 
plified versions of GR. As an example, some current core 
collapse codes (see e.g. @) use the ADM formulation of 
GR, along with the approximation of a conformally flat 
spatial metric [l^, [HI • Our splitting method provides a 
simple alternative way to add GR, without making any 
approximations. It should also be possible to split off 
and evolve waves equations about a curved background 
(such as Schwarzschild or Kerr) , which would be useful in 
simulating accretion disks or white dwarfs and neutron 
stars orbiting a supermassive black hole. 



II. FIELD EQUATIONS 

We will briefly review the GH formulation of GR - see 
e-g- [I3j Ell more details. We start with the metric 



ds^ = gabdx^'dx 



(1) 



The Einstein equations Gat ~ ^T^Tab can be written in 
trace reversed form: 



Rab = 47r(2Tafc - gabT), 



(2) 



with the Ricci tensor Rab given in terms of the connection 
coefRcients rjjj,: 

Rab = Tafc.c ~ ^'cb.a + ^'ab^dc ~ ^cb^da- (3) 



The lowered-indices connection coefficients are given by 
derivatives of the metric: 



1 



:{gab,c + 9ac,b — gbc,a), 



(4) 



and the metric can be used to raise the first index or get 
a contracted form: 



ra ^ad-r^ 
be — 9 ^ dbc, 



9 i a6c- 



(5) 
(6) 
(7) 



The contracted connection coefficients are equivalent 
to minus a wave operator acting on the coordinates: 
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g'"=VbVcx" = ax" = -^^ 



(8) 
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Choquet-Bruhat showed that if one picks harmonic 
coordinates Dx" = 0, then the Ricci tensor can be ex- 
panded and the Einstein equations ([2]) transformed into: 

9"'9ab,cd + 25^^,a56)d,c + Zr^.F^, (9) 

= -87r(2rafc - 9abT). 

The highest order term in Eq. ([9|) is the curved space 
wave equation g'^'^gab,cd: so this places the Einstein equa- 
tions in a manifestly hyperbolic form. 

A wider range of options are available however - one 
has the freedom in General Relativity to pick any set of 
coordinates. As shown by Friedrich [15j| and Garfinkle 
[l6t this freedom can be alternatively expressed by pick- 
ing source functions H°'(x, t) to drive the wave equations 
for the coordinates: 



Dx" = H''{x,t) 



(10) 



Using these source functions we can decompose Eq. ([2]) 
in the standard GH formulation: 

g'^gab,cd + 2g'=\^agb)d,c + 2i?(a,fc) - 2H,t:, (ii) 

+2r^,r:^„ = ^8n{2Tab - gabT). 

We will now modify this equation in order to explicitly 
separate off the wave equations, and convert the rest of 
the Ricci tensor into sources. We will proceed generally 
at first by considering metric perturbations fab away from 
a general background metric gab, (such as a Minkowski 
or Schwarzschild spacetimc): 

gab = gab + fab- (12) 

We then refer to the raised indices perturbations by h°'^: 

ab , f^ab (^3) 



9 
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(note that the /i"'' perturbations are determined by the 
background metric gab and the lowered indices perturba- 
tions fab)- 

One can then split the second order piece of Eq. (fTTj) 
into: 

g'^gab.cd = if" + h'''){gab,cd + fab.cd)- (14) 

We want to separate out a scalar wave equation for each 
of the metric perturbations fab on the background geom- 
etry gab- 



g'^'^Vc^dfab = g'^'^fab,cd — fab,e^'^- 



cd 4 



(15) 



We thus need to subtract fab,e^'^ from both sides of Eq. 
(fTTj) in order to get D/ab on the LHS. We collect the rest 
of the terms on the LHS of Eq. (fTT|) into a source Sab 
which we will move to the RHS: 



Sab — —g'^ gab.cd — h'^ [gab.cd + fab,cd) (16) 

- 2H^a.b) + ^H.r:, 

"^^ cb^ da ^ Jab,e^ - 



o cd 

{,agb)d,, 



This allows us to rewrite the standard GH equations (fTTj) 
in the simple form: 



Ofab — Sab ~ 87r(2rafc — gabT). 



(17) 



In Eqs. (jl6|) and (|17p we give a splitting for metric per- 
turbations on a general background geometry (jab- Here 
we specialize to the case examined in this paper: we set 
the background to be flat and use Cartesian coordinates 
Tjab = diag(-l, 1, 1, 1): 



gab — Vab + fab- 



(18) 



In this case the coordinate source functions H°' are zero 
and the Sab source terms simplify to: 

Sab = -h'"^fab,cd - '2.g'"^{^agb)d,c ~ '^'^tb^da (19) 

(where all of the metric derivatives now stem from the 
fab functions). 



III. NUMERICAL SIMULATION 
A. Matter EOM 

In order to test our splitting (jlGlllTp of the Einstein 
equations and see if they allow for stable and accurate 
numerical evolutions, we build a simplistic model of a bi- 
nary neutron star inspiral and merger. We construct and 
evolve a stress energy tensor Tab by fiat during the sim- 
ulation, which sources the gravitational fields fab (which 
in turn give rise to the corrective sources Sab)- Note that 
we are not using T°-^.b = to generate realistic equations 
of motion for the matter, as the main goal for this (sin- 
gle processor) code is to test the response of the gravita- 
tional fields fab to the rapid accelerations of dense matter 
sources. Later versions could use quas iequilibrium binary 
neutron star initial data as in 13, IS | , or fully evolve the 
matter equations of motion as in 19 1. 

The neutron stars are constructed as rigid polytropes 
with a compactness in the range of: M/R ~ 0.1 — 0.3. We 
choose to drive these density profiles on a quasi-circular 
inspiral path as found by Peters and Mathews [13, [2l[ , 
which captures the leading order radiation reaction ef- 
fects. The inspiral is parameterized by: 



a{t) flo 1 



td: 



■ecay 



1/4 



(20) 



where a{t) is the semimajor axis as a function of time, 
flo is the initial semimajor axis (between the centers of 
masses of the bodies), and the decay time is given by: 



td 



at 



ecay 



64 Af3 ' 



(21) 



where M is the total mass of the (equal mass) neutron 
stars. The instantaneous coordinate velocities of the 
stars are simply given by the Keplerian velocities of a 
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binary with the same mass and separation. Note that 
eventually the separation a{t) will decrease enough that 
the stars begin to overlap: here the stress energy tensor 
is determined by a simple superposition of the individ- 
ual stars densities, with the separation quickly falling to 
zero and the velocities also modulated to zero. This is, 
of course, not a particularly realistic model of a binary 
neutron star merger, the point is just to demonstrate 
that the gravitational fields fab still evolve stably in this 
strong field and highly dynamical setting. 

Having determined the bulk motion of the stars, we 
need the distributions of their fluid variables in order to 
build the stress energy tensor. We choose the standard 
perfect fluid form: 



T 



ab 



{p(l+e)+p)u°'v!' +pg 



ab 



(22) 



with rest mass density p, specific internal energy e, pres- 
sure p, and four velocities The internal energy and 
pressure are given in terms of the density through a poly- 
tropic equation of state: 



K r-i 



(23) 
(24) 



where we choose F = 2, and a value for k such that the 
internal energy density is about 5% of p in the center of 
the star. 

The density p is in turn derived from a conserved 
"baryonic" density p* (see e.g. [22]): 



(25) 



This density has a number of convenient properties, in- 
cluding a flat space conservation law: 

dtp* = 9,(«V*), (26) 
which ensures that the total integrated baryonic mass: 



M* 



(27) 



is constant throughout the evolution. Wc thus choose 
initial radial profiles for p* for each star and hold these 
constant during the evolution. 



B. Field Evolution 

Having determined the form of the matter stress energy 
tensor, we now concentrate on how to evolve the gravi- 
tational fields fab in response via (fT7|) (which has been 
simplified by Eqs. (jlSlllQp ). We proceed in a slightly 
nonstandard way: given previous successes we choose to 
finite difference the scalar wave equation Ofab using im- 
plicit methods, while we use the standard explicit meth- 
ods to construct the source Sab (note that this is a partic- 
ular choice, but any other stable scheme to evolve wave 
equations should work as well). 



First consider the implicit finite differencing of Ofab- 
We use operator splitting to divide the 3-1-1 wave equa- 
tion into three 1+1 wave equations (see e.g. [l!])- For 
conciseness replace fab with ip and the entire source 
Sab — SiTT(2Tab — QabT) with T. Adopting an operator 
splitting strategy we divide: 



- + dl^P H 
into a set of 1-1-1 problems: 



1/3t, 



+ dl^j = 1/3t, 
-d^ + dli: = 1/3t, 



(28) 

(29) 
(30) 
(31) 



which are sequentially finite-differenced implicitly. We 
use a variation of the Crank-Nicholson method, modified 
for second time derivatives, to evolve these - for instance 
Eq. (gni) becomes: 



1 

Ai2 



(32) 



where the field is discretized at time step N and spatial 
mesh location i. We have also included generalized coeffi- 
cients cjv+i and cm~i (with cm+i + cn-i — 1) so that we 
can transition from Crank-Nicholson CAT+i — cn~i = 1/2 
to a fully implicit method cjv+i = 1, cn-i — 0, or some- 
where in between. More implicit splitting allows us to 
dissipate high frequency noise if needed. Boundary con- 
ditions are added to Eq. which is then solved by 
inverting the resulting tridiagonal matrix. We note that 
in general implicit methods allow one to take as large a 
time step At as desired and still maintain stability, how- 
ever for reasons of accuracy we take time steps that are 
about the same size an explicit scheme takes in order to 
satisfy the Courant stability condition. 

This operator-splitting implicit method is combined 
with a Fixed Mesh Refinement (FMR) grid structure so 
that a high resolution mesh is available near the origin 
to resolve the stars, while also extending through suc- 
cessive lower resolution meshes out into the wave zone. 
Sommerfeld outgoing wave boundary conditions are ap- 
plied to the largest and coarsest mesh, which is updated 
first. Successively finer meshes are then updated, with 
their boundary conditions found via interpolation from 
the next coarsest mesh. After all the meshes are updated 
one time step the coarser mesh points are replaced with 
the finer mesh's field values wherever there is overlap (all 
meshes are updated each time step - no subcycling in 
time is used). 

This leaves the evaluation of the metric source Sab 
given by Eq. (jl9p . The spatial derivatives are computed 
in the standard second order explicit way, but evaluating 
the time derivatives is more subtle as the updated field 
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FIG. 1: The I = m = 2 mode of ^'4 as a function of time, 
demonstrating the standard chirp waveform. 

values /^^^ will not be known until Dfab is calculated. 
We choose to solve this iteratively: for the first pass the 
first and second time derivatives contained in (|19p are 
evaluated using the current and two previous time steps 
/afci fab^^^ fab "^ ^ this allows for a prehminary value of 
/j^'*'^ to be found. The second iterative pass then uses 
the preliminary /^^^ to calculate centered time deriva- 
tives in (jl9p which are used to complete the time step. 

Finally initial data is needed - we choose very simple 
initial data: fab{t = 0) = 0. There is thus a large amount 
of noise initially as the fields respond to the matter and 
then settle down to their correct physical solutions, with 
the transients propagating off the grid within a fraction of 
an orbit. For values of neutron star compactness much 
larger than M/R ~ 0.1 it is useful to also add a large 
amount of dissipation initially, and then linearly decrease 
the dissipation to a small or zero value over a fraction of 
an orbit. Otherwise the large amplitudes of the initial 
transient spikes can lead to infinities in the Sab sources 
(for example if |/oo| > 1), thus crashing the code. Using 
physical initial data could also solve this, but it is encour- 
aging that the code still rapidly converges to the correct 
solution when given initial data which completely ignores 
the self gravity of the stars (as has been seen elsewhere, 
see e.g. [24|). 

C. Results 

We find that our method of splitting off flat space 
wave equations for metric perturbations and turning the 
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rest of the Ricci tensor into sources provides an effec- 
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